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We study interacting bosons in a two dimensional square bipartite optical lattice. By focusing 
on the regime where the first three excited bands are nearly degenerate we derive a three orbital 
tight-binding model which captures the most relevant features of the bandstructure when the first 
excited p-bands in another sublattice are nearly degenerate with s-band of the other sublattice. In 
addition, we also derive a corresponding generalized Bose-Hubbard model and solve it numerically 
under different situations, both with and without a confining trap. It is especially found that the 
hybridization between sublattices can strongly influence the phase diagrams and in a trap enable 
even appearances of condensed phases intersecting the same Mott insulating plateaus. 
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I. INTRODUCTION 

The understanding that Hubbard models can be re- 
alized with ultracold atoms in optical lattices |l| has 
stimulated extensive effort to explore different aspects 
of quantum many-body physics in optical lattices [2|, |3| . 
The early works focused on the lowest energy band and 
in a pioneering experiment by Greiner et al. [4| the Mott- 
superfluid transition with ultracold bosons was observed. 
More recently, experimental groups have started to probe 
the properties of ultracold atoms under circumstances 
where the excited energy bands @ can no longer be ig- 
nored. This is most relevant since it has been demon- 
strated that the emerging multi-orbital effects can in- 
deed have crucial effects also on the ground state phase 
diagrams |6|. These excited bands can become impor- 
tant either when the atom-atom interactions become very 
large |7H14|. or when atoms are deliberately prepared 
on the excited bands. Such 'out of equilibrium' state 
preparation has been established by using accelerating 
lattices [15| or Raman transitions between bands |l6j . In 
the realm of these new experiments, one hopes to explore 
the regime where meta stable excited many-body states 
show very different properties from those of the ground 
state [12-^] ■ 

The experiment most closely relevant for our purposes 
is the one by Wirth et al. 2J]. Bosonic atoms were pre- 



pared in the ground state of a bipartite optical lattice and 
then the lattice was suddenly changed so that the initial 
ground state band atoms became (quasi) degenerate with 
a set of other bands which were initially separated by 
a large band gap. This process drove atoms into bands 
with non-trivial orbital properties and enabled the obser- 
vation of superfluidity on these so called p-bands. This 
experiment was followed by others |25l , |26| where uncon- 
ventional superfluidity was observed in the even more 
excited /-bands. 

Motivated by these experiments and especially on 
the aspects of the physics when different bands become 



degenerate we study multi-band bosons in a bipartite 
square lattice when bands cross. Such band cros sing can 
imply topologically non-trivial bandstructures [23, ^^. 
In principle, with the help of artificial gauge fields such 
bandstructures can also be engineered on the lowest 
band [23, [30| , but they might be easier to engineer in the 
excited bands were artificial gauge fields may become un- 
necessary. For example, in a square bipartite lattice the 
bandstructure can be composed of flat bands intersect- 
ing Dirac cones which, on the one hand, have interesting 
analogs with graphene physics, but the fiat bands also 
have novel influences on the dynamical properties of the 
gas [3l|, [33 • Physics of Dirac fermions have been studied 
in square optical lattices also in the absence of the fiat 
band 0]. 

As in the experiment by Wirth et al. ^if], we con- 
sider a bipartite square lattice of deep ^-sites and more 
shallow S-sites which however have a higher energy off- 
set. Under such circumstances, the excited (localized) 
states in A sites can become resonant with the ground 
states in B sites. When this happens, the p-bands can 
be strongly hybridized with the d-band. For vanishing 
atom-atom interaction, most of the relevant physics is 
captured by a tight-binding (TB) model, which predicts 
the existence of Dirac points and a flat band. Proceeding 
by adding atom-atom interactions we derive a generalized 
multi-band Bose-Hubbard model. We solve this theory 
from weak to strong interactions as well as in a trap. 
The calculated solutions reveal transitions from incom- 
pressible Mott insulators to condensed phases, but due to 
different atom-atom interactions the Mott lobes can be 
very dissimilar from those predicted by the usual single- 
band Bose-Hubbard model. Furthermore, the solution in 
a trap reveals the possibility that condensed states in dif- 
ferent sublattices occur in different regions of the trap. 
Our findings complement some other very recent ones, 
like Ref. [3^] where p-band bosons in a shallow bipar- 
tite optical lattice in terms of a nonlinear boson model 
is studied, and the work [35| analyzing the band struc- 



ture renormalized by the presence of interactions and the 
condensate in the broken symmetry phase. Finahy, Sun 
et al. [36 .;] also derived a fermionic tight-binding model 
which is quite similar to the one used by us. 

The paper is organized as follows. We begin by out- 
lining the theory relevant for our purposes in Sec. |IT1 In 
particular, Sec. Ill A] presents the tight-binding model to 
describe the ideal gas of atoms and in Sec. lIIBl we extend 
the model to include atom-atom interactions. In Sec. IIIII 
the generalized Bose-Hubbard model is solved within the 
Gutzwiller ansatz approach, and in Sec. IIII Al we discuss 
the solution in a harmonic trap. We end with a few con- 
cluding remarks in Sec. IIVI 



II. THEORETICAL FORMULATION 



A. Ideal system 




We will assume a two-dimensional lattice potential 
similar to the one used in the experiments by Wirth et 
al. (H; 
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where Vq is the lattice depth, k the lattice wave number, 
r] accounts for a small difference in the powers directed 
to different interferometer branches, e characterizes the 
power reduction in the retro-reflected beams due to im- 
perfect optics, and the angle a tunes the anisotropy intro- 
duced if e 7^ 1 . The angle 9 sets a relative phase between 
the two standing waves. £, y, and z are the unit vectors in 
the respective directions. Furthermore, the transverse z- 
direction has been reduced due to tight confinement. We 
will mostly consider a symmetric lattice with e — rj = 1, 
and cos (a) = e, but since different parameter choices can 
break the p-band degeneracies we allow for such possibil- 
ities as well. In Fig. [T] we show an example of a unit 
cell of this potential. Generally, the lattice is a bipar- 
tite square lattice where the two sublattices have lattice 
sites of different depths. Here we are interested in the 
parameter regime where the ground state in the shallow 
sites is quasi resonant with the first excited states of the 
deep sites. The resulting bandstructure of the regime we 
are interested in is depicted in Fig. [5] (a) and (b). Here, 
and in the following, we scale the energies in terms of 
the recoil energy Er = Ti (27r/A)^/2m of the atoms with 
mass m to absorb a photon of wavelength A. In partic- 
ular. Fig. [5] is calculated for Vq = IQEr. In this region, 
the two lowest excited p-bands become degenerate with 
the d-band. When this happens non-trivial bandstruc- 
tures with Dirac points emerge. Furthermore, one of the 
bands is almost flat suggesting that interactions play a 
larger role for atoms prepared in this band. 



FIG. 1. (Color online) The symmetric lattice potential with 
Vo ~ 10 Er over one unit cell. The parameters were chosen 
as e = 7j = 1, a = 0, and 9/tt = 0.556. xr and yR refer to 
coordinate axis rotated by 7r/4 with respect to the laboratory 
axes X and y. The shallow B-site is in the center while the 
deeper ^-sites are in the corners. Distance, A/2, between A- 
and B-sites was taken as a unit of length. 



Restricting our analysis to the three bands of Fig. [2l 
i.e. the localized ground state in the shallow S-sites and 
the first two excited states in deep ^-sites, we obtain an 
effective theory in terms of three different orbitals. In 
the absence of an external trap we can write the ideal 
gas Hamiltonian in momentum space as 



H = J2'PiHik)ct>i., 



where 



4^ 



(2) 



(3) 



describes the three types of orbitals included in our the- 
ory. There is an s-like orbital in the shallow i3-sites, ipf]^ 
and p-like x- and y-orbitals in the deep .A-sites, i/'^k and 

V'^k respectively. When the energy of the s-orbital in the 
S-sites is close to the energy of the p-orbitals in the A- 
sites, the dominant tunneling process is the one hybridiz- 
ing orbitals in different sublattices. This involves nearest 
neighbors and lower barrier height for tunneling while 
other tunneling processes require couplings over larger 
distances and are therefore greatly suppressed. Thus, for 
sufficiently deep lattices we can ignore tunnelings within 
A- or ;B-sites. On the other hand, since they only involve 
single particle physics, our theory can naturally include 
next nearest neighbor tunnelings easily when those are 
required. 



In momentum space, this results in a TB model 
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(4) 

whose parameters can be deduced from the exact band 
structure calculations, see Fig. [2l In the next section, 
this model will also be given in position space. One con- 
sequence of the hybridization can be seen in how the 
orbital character of the system enters for example in the 
sin-terms in the above TB model. Hopping occurs be- 
tween s- and p-orbitals, which implies that the tunneling 
coefficient alternates signs between neighboring sites giv- 
ing rise to a sin- rather than a cos-dispersion. In order 
to simplify notations, we choose our zero energy level 
to be the energy of the s-orbital in the K-sites. Since 
only nearest neighbor tunneling processes are included, 
the momentum dependence disappears from the diagonal 
terms and we have Sf (k) = 0, -B^(k) = E^ = 5/2, and 
E^{\.) = E^ ^ -5/2. 

We note that a somewhat related TB model was also 
derived by Sun et al. [3a|. However, in that model the un- 
derlying lattice potential was different and the p-orbitals 
were degenerate while in our case they can be different 
to account for the possible anisotropy of ^-sites. This 
anisotropy was indeed an important ingredient in the ex- 
periment by Wirth et al. [2J] . In the symmetric case with 
(5 = the lowest energy state of the TB model is 4-fold 
degenerate, but this degeneracy is lifted as soon as 5 ^ 
so that the minima is only two-fold degenerate. 

Furthermore, in the symmetric case with (5 = the sin- 
dispersions give rise to Dirac points at the origin as well 
as on the edges of the first Brillouin zone at (±7r/-\/2, 0) 
and (0, ±7r/v2)- A non-zero detuning 5 implies an effec- 
tive mass term that split the Dirac point degeneracies. 
Similarly, in graphene the relativistic electrons become 
massive when the symmetry between the corresponding 
two triangular sublattices is broken [23] ■ Contrary to 
graphene, rather than having a two-level structure, the 
present model has three bands and the Berry phase as 
a Dirac point is encircled vanishes. The additional level 
appears as a flat band sandwiched between the other two 
bands. 

H{k) of Eq. (U]) has the same structure as the Hamilto- 
nian for a A-scheme frequently occurring in light-matter 
interaction models in quantum optics, and we can di- 
rectly conclude that states of the flat band correspond to 
dark states with zero energy. These eigenstates are super- 
positions of p-orbitals and have a vanishing amplitude of 
being in the ( "excited" ) s-state in S-sites [33| ■ With this 
in mind, by considering anisotropic lattices (i^^ ^ ^yy) 
we notice that it would be possible to apply various exam- 
ples of complete or fractional stimulated Raman adiabatic 
passage schemes [38^ to prepare specific orbital states for 
the atoms. Intriguingly, the Hamiltonian in Eq. Q also 
has a clear connection to spin-orbit coupled systems. In 
the long wavelength limit we can expand the trigonomet- 



ric functions and flnd that the coupling between orbitals 
is linearly proportional to momentum [39l44l| . Usually 
spin-orbit coupling in ultracold atom systems is gener- 
ated between different atomic hyperflne states [411, '43] . 
Here the internal states of the atoms are not effected, 
but the spin-orbit-like coupling is a bandstructure effect 
that occurs between different orbitals. 

In Fig. [21(c) and (d) we demonstrate that the TB model 
above is indeed a good approximation close to band de- 
generacy by comparing it with the numerically calculated 
bandstructure, plots (a) and (b). As can be seen, for the 
symmetric lattice it reproduces the main features of the 
real bandstructure very well. Corrections beyond nearest 
neighbor hopping terms is seen to give rise to higher or- 
der variations in the dispersions mostly clear in the flat 
band. The tunneling coefficients t^^ and t:^^ have been 
extracted from the band widths of the numerically ob- 
tained bands. While our model does works well close to 
resonance, it should be kept in mind that generally the 
real bandstructure is more complicated and more tunnel- 
ing processes might have to be included in the theory. 



B. Interacting system 

In the previous section the ideal gas theory was de- 
rived and we now proceed by adding the atom-atom in- 
teractions. For ultracold atoms, interactions can be well 
modeled by contact interactions, 



C/= I / dri/it(r)V?(r)^(r)v3(r) 



(5) 



In a deep lattice the field operator ijj{r) is naturally 
expanded in terms of the localized orbitals described 
by the Wannier wave- functions w^{x,y), w:^{x,y), and 
wf{x,y). That is, we truncate the Hilbert space to con- 
tain only the three most relevant bands, i.e. the ex- 
pansion is restricted to ,B-sites' s-orbitals and ^-sites' 
p-orbitals. 

In the usual way, we limit the interaction to include 
only the dominant onsite terms. The strengths of vari- 
ous interactions are proportional to the scattering length, 
but their relative magnitudes depend on the orbital 
wavefunctions. To estimate these strengths we approxi- 
mate the onsite orbitals with harmonic oscillator wave- 
functions and in this way can analytically solve the inte- 
grals describing interaction between cc-orbitals in ^-sites 



Uxx = ^0 / dxdy\w^{x,y)\'^, 
between y-orbitals in ^-sites 

Uyy = f/o / dxdy\Wy{x,y)f, 
between x- and y-orbitals in ^-sites 

Uxy = Uo dxdy\w^{x,y)\'^\Wy{x,y)\^, 



(6) 



(7) 



(8) 




FIG. 2. (Color online) Dispersions of the three lowest excited bands, (a) and (b) are numerically calculated for a lattice with 
Vo — 10 Er, e = 1] = 1, a = 0, and 9/tt = 0.556. (a) shows the lowest excited band while (b) shows all three excited bands in 
the same plot, (c) and (d) are calculated with the TB model with parameters t^^ — tyy = 0.06485i5H and 5 — Q. 



and finally between s-orbitals in S-sites 



UsB = Ua dxdy\wf{x,y)\' 



(9) 



We take that the remaining prefactor C/q is tunable either 
by changing the lattice depth or by changing the effec- 
tive scattering length. In the harmonic approximation 
Uxy = Uxx/'i- This condition can sometimes lead to ac- 
cidental degeneracies, which are removed as soon as the 
condition is broken [2J]- However, in this work this does 
not play a major role. Since the shallow sites are wider 
than the deep sites, their orbitals are also more extended. 
This implies that Usb is often surprisingly close to the 
values of Uxx and Uyy even though these involve wider 
excited state orbitals. For concreteness, in the following 
we choose the lattice depth as Vq = 10 Ej^ in which case 
it turns out that Uxx = Uyy ~ 0.95 Usb- 

With the above introduced interaction strengths, we 
are now in a position to write down a many-body Hamil- 
tonian describing multi-orbital bosons in a bipartite op- 
tical lattice. The corresponding Hamiltonian takes the 
form 



where 
Ho 
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(11) 



Ht — Hq + Hj^js + Hi,j,, 



(10) 



describes the energy offsets and nearest-neighbor tunnel- 
ing giving rise to hybridization between orbitals. Here 
i = {ix,iy) labels the lattice sites and /i is the chemical 
potential, n■^l^ fiyi, and n^ ; are the number operators for 
X- and y-orbitals in an y^-site i and s-orbitals in a ,B-site 
i. The notation j^^ (j^_ ) indicates a nearest neighbor of 
i — {ixjiy) to the right (left) in the direction /3 e {x,y}. 
For example, j^,^ = {ix + l,iy) while jx = [ix - l,«i;)- 
Finally, h.c. indicates the hermitian conjugate. The hop- 
ping term must be written in this way since in this case 
tunneling is sensitive to the left and right difference. In- 



tuitively this is easy to understand by considering a p- 
orbital with a node. This orbital wave function changes 
sign as one moves along to axis towards the neighboring 
site with s-orbital wavefunction. The overlap of these two 
wavefunctions is predominantly positive if the neighbor 
is to the left (for example) , but predominantly negative if 
it is to the right. Note how such "space-dependence" in 
the hopping term also appears in lattice models exposed 
to (synthetic) magnetic fields [4l| . The tunneling param- 
eters t-^p denotes the strength of tunneling of a-orbitals 
in the A sublattice in the direction j3 into the nearest 
neighbor s-orbital in the B sublattice. In the theory used 
here t-^y = t^ = 0. In the momentum representation, 
the term Hq corresponds to the TB Hamiltonian encoun- 
tered in the previous subsection. 

The remaining terms describe interactions. 






1 



(12) 
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accounts for the interactions in the S-sites and 
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(13) 
for the interactions within the ^-sites. This term is some- 
what more complicated than the corresponding term in 
the shallow sites since x- and y-orbitals interact and (for 
bosons) can change into one another. Finally, we note 
that when (5 = the total Hamiltonian supports a sym- 
metry corresponding to swapping of x- and y- flavor atoms 
in the .A-sites. 



III. GUTZWILLER RESULTS 

The Gutzwiller ansatz [43| for the many-body wave 
functions provides a reasonably accurate description of 
interacting bosonic systems, especially in dimensions 
D > 1. Due to the bipartite lattice and multiple flavors in 
one sublattice, our case is somewhat more complex than 
the usual Bose-Hubbard model. The Gutzwiller ansatz 
we use is given by 



HE 
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(14) 
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The expansion coefficients a K — a K 



di) 



and b L are the 



Gutzwiller amplitudes of the corresponding on-site Fock 
state. For our purposes, in the yi-sites the relevant sub- 
space is spanned by the Fock states of the form |n-^) — 
|n;^,n^), where n^ is the occupation number of the a- 
orbital. In the S-sites the wave function is expanded in 



terms of Fock states |nf) associated with s-orbitals. The 
Gutzwiller ansatz captures the onsite physics exactly, but 
ignores some correlations between sites. In the limit 
where the onsite wave function are taken to be coher- 
ent states it recovers the Gross-Pitaevskii limit of weakly 
interacting bosons. This limit is approached as interac- 
tions relative to kinetic energy become small. In the limit 
of strong interactions, the Gutzwiller ansatz can predict 
different insulating phases. Depending on the problem, 
the insulating states predicted by the Gutzwiller ansatz 
can be degenerate and these degeneracies can in princi- 
ple be broken due to the weak inter-site correlations not 
encountered for in this approach. This was demonstrated 
for the square and cubic lattices by treating the kinetic 
energy terms as perturbations [23| . 

Calculating the energy expectation value ('0|iJT|V')i 
using the ansatz in Eq. p4p . gives us an energy functional 

in terms of the unknown (complex) amplitudes a j^ and 
h B . This energy functional must then be minimized to 
flnd the ground state. Even though this functional is 
very complex and the minimization is not always easy, 
we have found that standard conjugate gradient meth- 
ods work with few caveats. First, the energy functional 
can have many local minima into which the minimiza- 
tion algorithm can become stuck and consequently fail 
to converge into the global minimum. In order to build 
up confidence in the results it is important to try different 
initial states. Second, the minimization algorithm might 
have trouble in converging to the correct phase ordering. 
For example, complex amplitudes give rise to different 
phase factors in the condensate order parameters and in 
an energy minima these phase factors should be prop- 
erly ordered throughout the lattice [21]. If the conjugate 
gradient method is used as a black box, it might not con- 
verge to optimal phase ordering. To get around this, it 
is important to impose different orderings into the ini- 
tial state of the minimization routines and finally pick 
the solution that has the lowest energy. In the absence 
of a trap we find the solution in a 4 x 4 lattice where 
each sublattice has 8-sites. We use periodic boundary 
conditions and choose to truncate the Fock state expan- 
sion of the Gutzwiller ansatz so that the maximum onsite 
occupation number is 8. 

In the superfiuid region wc flnd that the ground state 
phases of the condensate order parameters are arranged 
in the same way as discussed by Wirth et al. [2J] for 
an isotropic lattice. Here, the phase of the condensate 
order parameter in the S-sites changes by ±27r as one 
moves around S-sublattice plaquettes. Neighboring pla- 
quettes have an opposite phase winding. In the yl-sites, 
the onsite order parameters are superpositions of x- and 
2/-orbitals. These superposition are vortex- like states pro- 
portional to e"^ {x ± iy) and the vorticity has an opposite 
sign in neighboring y^-sites so that onsite angular mo- 
menta are ordered "anti-ferromagnetically" . Similarly to 
the ,B-sites, the phase of the phase factor e**^ varies by 
±27r as one travels around ^-site plaquettes and neigh- 



boring plaquettes have an opposite winding of this phase 
factor. Far in the superfluid phase where the onsite states 
can be approximated by coherent states, atoms in the 
y^-sites can be pictured as clockwise or anti-clockwise ro- 
tating condensates with a quantization (i^) = ±1 where 
Lz is the angular momentum operator in the transverse 
z-direction. 

Note that the swapping symmetry of x- and y-flavor 
atoms implies a flip of the vorticity in each site. Closer to 
the insulating phases where interaction begins to domi- 
nate, the picture is more complex and the onsite x- and 
y-flavor atoms can become highly entangled. While the 
Gutzwiller ansatz ([T4| is not able to predict inter-site 
entanglement it indeed captures such intra-site entangle- 
ment. As an example, looking at the Mott insulating 
phase with n-^ = 3 atoms in the .A-sites, the Gutzwiller 
method gives a degenerate ground states in the yl-sites. 
For example, the states with \ili)L — Y\iizj,a'^Q\'i,^)i + 



''1,2 
*3,0 



l,2)i or \i^)n = ni6^<3lO,3)i 



^0,3 



0,31 

(i) 
2,1 



,« 



2,l)i, with 



AjQ = 0^3 « 0.6 and a\2 = ^21^ —0.8 are degenerate. 
As discussed above, in the Gutzwiller method these two 
states are decoupled in the insulating phase and breaking 
the degeneracies requires improved ansatz and/or higher 
order perturbation theory in tunneling. It is clear that 
these two examples of insulating states are not eigen- 
states of Lz- 

We show an example of the magnitudes of the rele- 
vant observables in the phase diagram for the isotropic 
case with degenerate p-orbitals in Fig. [31 As is clear, 
the phase diagram is very different from the usual se- 
quence of ever lower Mott-lobes corresponding to higher 
onsite atom numbers [44j . In our case there are insu- 



lating states with integer occupation numbers, but since 
interactions in different sublattices are different and the 
other sublattice has several flavors the positions of the 
boundaries for different Mott-states are not expected to 
be in same positions for different sublattices in the limit 
of weak tunneling. The hybridization of orbitals in dif- 
ferent sublattices complicates the picture further. 

This interplay between sublattices gives rise to super- 
fluid "fingers" extending into the region where each sub- 
lattice alone would be expected to be in a Mott insulator. 
For example, yB-sites make a transition from 1 atom per 
site to 2 atoms per site at /i/C/gg = 1. This is apparent in 
the order parameter (V'f ;)^ being non-zero in the narrow 
region around ^jl/Usb = 1 even when tunneling becomes 
weak. With these parameters and weak tunneling the 
.4-sites are expected to be in an insulating state with 2 
atoms per site (|n^ — l,n^ — 1)), but coupling with the 
condensate order parameter in the ;B-sites can induce a 
non-zero order parameters {ip0i)- Similar observations 
apply around ii/Uss ~ 1-25 where the ^-sites undergo a 
transition to 3 atoms per site. This transition can induce 
a non-zero condensate order parameter in the S-sites. 

It should be noted that the number fluctuations in x- 
and j/-flavors in the ^-sites can be non-zero even in Mott 
insulating regions. For example, the Mott insulating 



state with 3 atoms in the ^-sites is a superposition of dif- 
ferent basis states with the total of 3 atom per site. Only 
the total number of atoms is fixed to an integer value. 
The local order parameter breaks the time-reversal sym- 
metry and the angular momentum in ^-sites is non-zero 
and equal to ±1 in the condensed region. The angular 
momentum in neighboring ^-sites point in opposite di- 
rections. The non-zero value of angular momentum in 
the condensed phase is not surprising since the interac- 
tion energy is minimized for onsite states with x±iy type 
vortex superpositions of p-orbitals [23, \^ . 

In Fig. m we show an example of the phase diagram for 
the anisotropic case with a p-orbital splitting 5/UsB = 1- 
In the superfluid regions both p-orbitals are non-zero, 
but the order parameter (and density) for the y-orbital 
is smaller in magnitude. In the Mott insulating regions 
with 1 or 2 atoms per site, the onsite interactions (with 
these parameters) are not strong enough to induce large 
fraction of atoms into the higher y-orbitals and therefore 
only the x-orbitals are substantially populated. However, 
as the atom number in the ^-sites increases to 3 or more 
also the y-orbital population becomes substantial. 

When we choose (5 7^ we break the degeneracy of 
the X- and j/-orbitals. In the limit of zero tunneling we 
expect that if splitting becomes in some sense large rel- 
ative to onsite interactions, atoms would prefer to reside 
on the s-orbital only. It is easy to show that with 2 
atoms per ^-site, the transition occurs at (5 = Uxx/^- It 
is important to keep in mind that for the case of non- 
zero tunneling, the situation becomes much more com- 
plex and the results may actually depend on the system 
size. With the Gutzwiller ansatz we find that in the su- 
perfluid regime (we typically had t^x/UsB ~ 0-2 . . . 0.5), 
the onsite angular momentum (which vanishes if only one 
orbital is occupied) per particle is smoothly reduced from 
its value ±1 at J = to zero. This is demonstrated in 
Fig. [5] Vanishing onsite angular momentum is reached 
when 5/t-^^ ~ 2 which corresponds fairly well to what 
would be predicted from the onsite results with a par- 
ticle number fixed to an integer value (remember that 

UsB ~ Uxx) 

S/tt^ = {5/Usb) X (UsB/tt^) « I X (C/.g/e^) . (15) 

As expected, the onsite angular momentum also drops 
faster for larger t^^ I^sB since this implies smaller onsite 
interaction strengths. 

If we replace the operators with complex numbers ^q 
to derive a Gross-Pitaevskii equations for each orbital, 
we find that for the onsite problem the effective chemical 
potential and thus also the density of ?/-orbitals vanish 
when 5/2 = fj, — UxyUx at which point the density of 
the x-orbital is related to the chemical potential through 
Ux — {pL + 8/2)/Uxx- This implies that in this limit the 
transition from states with orbital angular momentum to 
pure a;-orbital condensate happens at 5c — {Uxx~Uxy)nx, 
where Ux = IV'^P- 
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FIG. 3. (Color online) Condensate order parameters and onsite atom numbers parametrized by the chemical potential and 



hybridization tunneling t 



AB _ ,AB 



tyy when p-orbitals are degenerate. (However, in order to make the plot clearer we did 



add a very small anisotropy oi 5 — 10"* to break the degeneracy of states in j4-sites with only one atom per site.) The left hand 
plots (a),(c), and (e) display condensate densities \{ii's,i)\ and |{'!/'^,i)| (/? € {x,y}), while (b), (d), and (f) show atom flavor 
densities nf j and n^^^. The roughness that is visible especially for higher chemical potentials indicates the level of numerical 
uncertainties in these regions. (In the Mott insulating region with n"^ = 1 we choose n;^; = 1, but since interactions do not 
contribute here other choices are also possible.) 




FIG. 4. (Color online) Condensate order parameters and onsite atom numbers parametrized by the chemical potential and 
hybridization tunneling t — t'^^ — tyy for the anisotropic case with 5/UsB ~ 1- The plots to the left, (a), (c), and (e), 
show condensate densities |{^/'^i)P and |{V'^_i)P (/3 £ {a;,j/}) while the ones to the right, (b), (d), and (f), display atom flavor 
densities ris^i and n^^,. Small amount of scatter visible especially in (e), is indicative of numerical uncertainties. 




FIG. 5. Angular moraentum per particle in the ^-sites as a 
function of energy difference 5 between p-orbitals. We choose 
ixx /UsB ~ 0.2, ^/t^x = 1] and t^^ as the unit of energy. 
(The staircase structure at larger 5 is due to numerical limi- 
tations in finding the global energy minimum for larger onsite 
atom numbers with a finite basis set.) 



A. Trapped system 



Typical experiments would most likely involve the 
presence of a confining trapping potential and for this 
reason it is important to also discuss the behavior with 
inhomogeneous density distributions. Our predictions for 
the phase diagram in a homogeneous system suggest an 
interesting possibility in a trap. Usually the solution of 
the Bose-Hubbard model in a trap gives rise "a wedding 
cake" structure where Mott plateaus corresponding to 
different integer fillings are sandwiched between super- 
fluid regions |46| . 

If we were to apply a local density approximation to 
our system, we could think of the chemical potential 
as a local quantity ^ == fJ-center - Vtrap{ix,iy), where 



Vt 



trap\^x: ^y 



) would typically be a harmonic trap. Travers- 



ing from the center of the cloud to its edge would corre- 
spond to moving in the phase diagram from some high 
value of i-i/UsB towards zero. If the starting point is in 
the Mott insulating phase we could indeed have a wed- 
ding cake structure for each sublattice, but their Mott 
plateaus do not always coincide. Furthermore, we can 
have situations when a condensate order parameter ap- 
pears inside the same Mott plateau. We will next demon- 
strate that these simple observations are valid in a trap 
also beyond the local density approximation. 

We can do this within the theoretical framework used 
so far, but replacing the chemical potential /i with 
fJ'center — Vtrap{ix, iy) in the Hamiltonian in Eq. dill) and 



then solving the problem with the trapping potential 

Vtrap{ix,ly)=l[{ix-iN^ + l)/2f + iiy-iNy + l)/2f] 

(16) 
with Nx and Ny being the number of sites along x and 
y respectively. (The minima of the harmonic potential 
is shifted to {{N^ + l)/2, {Ny + l)/2) since we choose 
ia € {1 . . .Na}.) As an example, we choose an isotropic 
lattice with t/Uss = 0.015 and the chemical potential 
in the center iJ-center/Usrs — 1-5 so that in the center 
of the cloud we expect the y^-sites to be in an insulating 
state with three atoms per site. The trap coefflcient 7 we 
choose in such a way that ^center — Vtrapiix, iy) becomes 
negative at the edge of the lattice so that the density 
vanishes there. 

We demonstrate the resulting ground state of the 
trapped bosons in Fig. [S] The bosons arrange them- 
selves into the familiar wedding cake structure with Mott- 
insulating regions separated by superfluid-regions. Re- 
markably, as suggested by the results in the absence 
of trapping potential, since our system has two differ- 
ent sublattices with different onsite interactions, super- 
fluid "rings" can occur in different locations for dif- 
ferent orbitals. For example, closest to the center we 
have a region where the ^-sites are Mott-insulators with 
n"^ — n^ -|- n^ = 3 while the S-sites are insulating with 
n^ = 2. The transition to n-^ = 2 phase occurs via a su- 
perfluid phase in the y^-sites. However, in this region the 
S-sites are still very small. Also, there is a condensed 



phase between regions with n — 2 and 



,B _ 



1 while 



the condensate order parameters in y^-sites are negligi- 
ble. Consequently, the physics predicted by using the 
theory without the trapping potential can also persist in 
trapped systems. 

Recently, the trapped system of p-band bosons in a 
square lattice was analyzed and it was found that the 
density of different x- or j/-orbital atoms were elongated 
in one direction and the symmetry of the confining trap 
was broken [i^. The present system is different due 
to the hybridization of s- and p-orbitals, which implies 
that the condensate cloud preserves the symmetry for an 
isotropic trap. On the other hand, if one prepares the 
system so that the tunneling coefficients t-^^ and id^ 
are unequal in magnitude, similar anisotropics might be 
expected also here. 



IV. CONCLUSIONS 

In this paper we have derived a TB model to describe 
ultracold atoms in a bipartite optical lattice with three 
hybridized orbitals. We have also solved the resulting 
generalized Bose-Hubbard model and found strong mod- 
ifications to the Mott insulator superfluid phase diagram 
which is found in the simplest lowest band Bose-Hubbard 
model. Novel phenomena was also demonstrated for the 
confined system that includes a harmonic trap. From 
that solution we found that the unusual phase diagram 
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FIG. 6. (Color online) Condensate and flavor densities in a trap. The left hand plots (a) and (b) give the condensate densities 
Ki/ifi)!^ and KV'^JP while (c) and (d) show atom flavor densities nf ; and n;^;. We choose t/UsB ~ 0.015, ficenter/UsB ~ 1-5, 
and 7 in such a way that the density vanishes at the edge of the lattice. Since the lattice is isotropic the densities for the 
j/-orbital are the same as for the a;-orbital and are not plotted here. The axes give the lattice sites in the two laboratory 
directions. (Plotted quantities are only deflned in their respective sublattices. However, to make the flgure clearer we flUed in 
the relevant values also to the other sublattice by taking the average over the 4 neighboring sites.) 



of the multi-band Bose-Hubbard model can be reflected 
as possessing non-trivial wedding cake structure of Mott 
insulating regions for different sublattices. In particular, 
a non-zero condensate order parameter in one sublattice 
can coexist with a Mott plateau in another sublattice 
and also appear inside the same Mott plateau. Such ef- 
fects are observable since Mott insulating regions can be 
detected in-situ and atoms in optical lattices can be ma- 
nipulated even at a single site resolution |47H50| . Fur- 
thermore, since different sublattices have different atom- 
atom interactions the states with more than one atom 
per site would generally give rise to different mean-field 
shifts if transitions to other hyperfine states are consid- 
ered. This suggest a possibility of addressing different 
sublattices with microwave fields of different frequencies, 
for example. 



In this paper we have not addressed the dynamical 
behavior of bosons in a bipartite lattice. However, using 
the theoretical framework derived here that would be not 
only doable, but also interesting since in the experiments 
conducted so far bosons have been initially prepared in 
an excited state whose dynamical behavior is poorly un- 
derstood. 
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